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ABSTRACT 


A computer  code  is  presented  for  the  simulation  of  turbulent  swirling 
flow  through  closed-conduit  flow  meters.  The  code  is  reasonably  fast  and 
accurate.  The  response  of  a flow  meter  to  changes  in  its  operating 
environment  is  easily  assessed.  Through  the  use  of  computer  plotting  routine 
a complete  picture  of  the  flowfield  through  a meter  can  be  obtained.  An 
example  flowfield  for  a target  meter  is  presented. 


1 . INTRODUCTION 

A computer  code  currently  in  use  at  NBS  has  proven  highly  useful  in 
determining  the  characteristics  of  closed-conduit  flow  meters  operating 
under  realistic  conditions  [1].  This  code,  called  TEACH,  has  been  developed 
at  the  Department  of  Mechanical  Engineering,  Imperial  College  of  Science 
and  Technology,  London,  England.  It  is  in  wide  use  throughout  the  United 
States  to  model  such  turbulent  flows  as  those  found  in  combustors  and  in 
nuclear  reactors  [2,  3].  The  TEACH  code  has  generally  been  found  to  give 
reasonably  accurate  results  within  reasonable  computing  times.  It  was  thus 
the  logical  candidate  to  pick  in  testing  the  feasibility  of  using  turbulence 
modeling  to  study  the  flow  through  various  types  of  closed-conduit  flow 
meters.  The  results  have  demonstrated  this  feasibility,  and  computer  simulation 
is  now  regarded  as  a vital  tool  in  NBS'  effort  to  improve  flow  metering 
practice . 

The  TEACH  code  can  be  used  to  obtain  steady  solutions  for  axisymmetric 
flows  with  or  without  swirl,  a swirl  subroutine  having  been  added  at  NBS.  The 
ke-model  is  employed  to  determine  the  turbulent  shear  stresses  [4].  The  resulting 
solution  is  in  terms  of  velocities,  pressure,  and  turbulence  intensity.  These 
can  be  computer  plotted  in  such  a way  as  to  give  a detailed  view  of  the  flowfield 
through  a metering  device.  Changes  in  this  flowfield  due  to  variations  in  such 
parameters  as  upstream  velocity  profile,  Reynolds  number,  or  swirl  number  are 
easily  assessed.  A target  meter  flowfield  is  presented  as  an  example. 

2.  BASIC  EQUATIONS 

The  momentum  and  continuity  equations  for  turbulent  incompressible 
stationary  flow  employed  in  the  ker-model  are 

p(q  * V)q  = - VP  - Vx[peff (Vxq) ] (1) 


V • q = 0 


(2) 


in  which  q = (U,  V,  W) , where  U,  V,  and  W are  the  mean  velocity  components  in 

the  radial,  axial,  and  tangential  directions,  respectively,  in  an  axisymmetric 

cylindrical  reference  frame;  P is  the  mean  pressure;  p is  the  constant  density; 

and  p cc  is  the  sum  of  the  laminar  and  turbulent  viscosities.  The  turbulent 
etr 

viscosity,  p , is  assumed  to  be  a function  of  the  turbulent  kinetic  energy, 
k,  and  dissipation  rate,  e,  of  k;  i.e., 


yt  = Cyp  k2/e  (3) 

where  the  empirical  constant  C is  0.09.  The  values  of  k and  e are  determined 

y 

by  the  following  transport  equations  written  in  axisymmetric  cylindrical 
coordinates  (r,  z)  : 
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where  the  empirical  constants  are  a,  = 1.0,  o = 1.3,  C = 1.0,  C = 1.44,  and 
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r 1 no.  , r _ 0 r ^ /U . 2 , ,3V.  2,  , r 9V  3U.2  .9W.2  .9  W.  . 2 
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I dr  r 9z  dr  dz  dz  dr  r 

The  origin  of  the  values  for  the  empirical  constants  is  discussed  in  [4]. 

3.  NUMERICAL  SCHEME 


The  TEACH  code  employs  the  SIMPLE  finite  difference  scheme.  This  scheme  is 
described  in  [5,  6],  while  the  code  itself  is  discussed  in  [2,  7].  Only  a brief 
review  of  TEACH  will  thus  be  given  here. 

The  TEACH  code  employs  a variably-spaced,  rectangular  grid.  The  radial  cell 
dimension  is  only  a function  of  radial  position,  and  the  axial  cell  dimension 
depends  only  on  axial  position.  Normal  velocities  are  specified  along  cell 
boundaries,  and  pressure,  tangential  velocity,  turbulent  kinetic  energy,  and 
dissipation  rate  are  defined  at  cell  centers.  The  system  of  equations 
(1)  - (5)  is  finite  differenced  and  solved  iteratively  by  underrelaxation.  The 
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nonlinear  convective  terms  in  these  equations  are  differenced  using  a hybrid 
scheme  comprised  of  central  and  conservative  upwind  differencing  [8].  At 
each  iterative  step,  the  finite  difference  analogs  of  equation  (1)  are  solved 
implicitly  using  a line-by-line  procedure  involving  a tridiagonal  matrix  solver. 

The  velocity  field  resulting  from  this  solution  of  the  momentum  equations  for  a 
given  iteration  will  probably  not  satisfy  the  continuity  equation  (2).  Using 
this  known  (but  probably  incorrect)  velocity  field,  equation  (2)  can  be  transformed 
into  an  equation  for  the  pressure  correction  required  to  make  the  velocity  field 
solenoidal.  This  pressure  correction  is  calculated  in  the  same  implicit 
manner  and  is  used  to  determine  the  final  velocity  field  for  a given  iteration. 
These  velocities  are  then  used  in  the  solution  (by  the  usual  implicit  scheme) 
of  equations  (4)  and  (5)  for  k and  e.  Finally,  equation  (3)  for  y^_  is  updated 
and  the  next  iteration  begins.  The  iterations  proceed  until  the  sum 
of  the  absolute  values  of  the  residual  mass  and  momentum  sources  is  less  than 
some  preset  small  value.  This  summation  is  performed  separately  for  the  mass 
flow  and  each  component  of  momentum,  with  the  resulting  sums  to  each  be  less 
than  some  small  percentage  of  the  mass  flow  and  largest  component  of 
momentum,  respectively,  at  the  inlet  to  the  computed  region. 

Boundary  conditions  must  be  specified  at  the  inlet  to  the  region  of 
flow  containing  the  metering  device.  The  required  conditions  are  radial 
distributions  of  U,  V,  W,  k,  and  e.  The  geometry  of  the  flow  meter  itself 
is  handled  by  not  allowing  fluid  to  cross  the  appropriate  grid  lines.  The 
condition  at  the  exit  of  the  computed  region  is  that  axial  gradients  are  zero. 

Shear  stresses  along  walls  are  based  on  the  logarithmic  law;  i.e.,  wall- 
function  method  of  [4].  Axisymmetry  is  enforced  at  the  centerline  (r  = 0)  of 
the  computed  section. 
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4.  COMPUTER  PLOTTING 


The  plotting  program  makes  extensive  use  of  four  CalComp  (California 
Computer  Products)  subroutines.  Three  of  these,  PLOT,  SYMBOL,  and  NUMBER, 
are  part  of  CalComp 's  basic  software  package  that  adapts  program  output 
for  plotting.  The  fourth  subroutine,  FLINE,  is  classified  as  a "functional" 
subroutine,  in  which  basic  subroutines  are  applied  to  more  specified  needs. 

Both  basic  and  functional  software  is  available  from  CalComp. 

The  PLOT  subroutine  converts  the  output  data  from  TEACH  (U,  V,  W,  P,  k,  e) 
into  hardware  codes  for  a pen  plotter.  SYMBOL  enables  printing  of  alphanumeric 
data,  such  as  titles  for  plots.  This  subroutine  also  creates  the  symbols  used 
to  plot  the  data  points.  NUMBER  is  used  to  provide  axis  value  labels  on  plots. 
The  FLINE  subroutine  has  the  capability  of  plotting  an  entire  array  of  data 
points  with  symbols,  and  then  drawing  straight  lines  or  smooth  curves  between 
successive  data  points.  A parabolic  curve-fit  is  generally  used  at  NBS . 

The  plotting  program  generally  runs  in  less  than  three  minutes  and  uses 
about  40K  words  of  storage  on  the  NBS  UNIVAC  1108.  About  1.5  hours  is  usually 
required  to  make  the  graphs  from  the  resulting  output.  When  a complete  set 
of  plots  is  generated,  approximately  9 meters  of  30  cm  paper  is  used.  Some 
examples  of  these  plots  for  a target  meter  are  presented  in  the  next  section. 

5.  EXAMPLE:  TARGET  METER 

As  an  example  of  what  can  be  done  with  TEACH  and  the  plotting  program, 
a target  meter  flowfield  is  presented  in  Fig.  1.  The  target  meter  is  merely 
a disc  of  zero  thickness  (T  = 0)  suspended  in  the  center  of  a circular  pipe  at  a 
distance  of  four  pipe  radii  from  the  inlet.  The  beta  ratio  (ratio  of  disc  radius 
to  pipe  radius)  is  0.5,  and  the  Reynolds  number  based  on  pipe  radius  and  mean 
inlet  velocity  is  2r10“*.  The  inlet  conditions  are  based  on  fully-developed  pipe 
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flow  at  this  Reynolds  number,  with  the  axial  velocity  profile  conforming  to  a 
standard  l/9-power  law.  A solid-body  rotation  is  superimposed  on  these 
inlet  conditions  in  order  to  assess  the  effects  of  swirl  on  the  meter.  The 
inlet  swirl  number,  which  represents  the  ratio  of  axial  flux  of  angular 
momentum  to  axial  flux  of  axial  momentum,  is  0.10.  All  quantities  are 
nondimensional , with  the  characteristic  length  being  pipe  radius  and  the 
characteristic  velocity  being  mean  inlet  velocity.  There  are  100  variably- 
spaced  grid  points  axially  and  20  variably-spaced  ones  radially.  It  required 
about  400  iterations  to  reach  the  final  solution,  which  took  about  40  minutes 
and  55K  words  of  storage  on  the  NBS  UNIVAC  1108. 

Figure  (la)  shows  pressure  and  axial  velocity  variations  through  the  meter. 
The  pressure  difference  across  the  disc  can  be  integrated  to  give  the  force 
on  the  disc,  which  is  related  to  flow  rate.  Figure  (lb)  gives  wall  shear 
stress  (x)  and  swirl  number  (S)  anywhere  in  the  meter.  Figure  (lc)  is  a 
streakline  (segmented  streamline)  plot,  with  each  line  segment  being  in  the 
local  flow  direction  at  a given  point.  The  recirculation  zone  along  the  pipe 
centerline  behind  the  disc  shows  clearly,  as  does  the  concentration  of  grid 
points  near  the  disc.  Figure  (Id)  presents  turbulence  intensity  based  on  local 
mean  velocity,  I,  along  the  pipe  wall  and  pipe  centerline.  The  flattened  peak 
along  the  centerline  results  from  cutting  off  a spike  which  occurs  at  the  end 
of  the  recirculation  zone  where  the  local  mean  velocity  goes  to  zero  (I  -*■  °°)  . 
Figure  (le)  shows  axial  velocity  profiles  at  several  critical  points  in  the  meter. 
Note  the  power-law  inlet  profile.  The  outlet  profile  has  not  yet  returned  to 
the  fully-developed  case.  Figure  (If)  presents  tangential  velocity  profiles  at 
the  same  locations  as  for  the  axial  ones.  The  inlet  profile  is  linear, 
denoting  solid-body  rotation. 
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6.  CONCLUSION 


The  feasibility  of  performing  computer  simulations  of  turbulent  swirling 
flow  through  closed-conduit  flow  meters  has  been  demonstrated.  A turbulence 
code,  such  as  the  TEACH  code  presented  here,  can,  along  with  proper  plotting 
routines,  define  entire  flowfields  quickly  and  cheaply.  It  is  thus  expected 
that  this  sort  of  simulation,  in  conjunction  with  well-chosen  validation 
experiments,  will  become  more  common  in  the  determination  of  the  performance 
of  flow  meters  and  associated  pipeline  elements,  such  as  upstream  flow  conditioners. 
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Fig.  1.  Target  meter  flowfield. 
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